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A new Self-Organized Criticality (Soc) model is introduced in the form of a Cellular Automaton (CA) for 
ion temperature gradient (ITG) mode driven turbulence in fusion plasmas. Main characteristics of the model 
are that it is constructed in terms of the actual physical variable, the ion temperature, and that the temporal 
evolution of the CA, which necessarily is in the form of rules, mimics actual physical processes as they are 
considered to be active in the system, i.e. a heating process and a local diffusive process that sets on if a 
threshold in the normalized ion temperature gradient R/Lt is exceeded. The model reaches the Soc state and 
yields ion temperature profiles of exponential shape, which exhibit very high stiffness, in that they basically 
are independent of the loading pattern applied. This implies that there is anomalous heat transport present 
in the system, despite the fact that diffusion at the local level is imposed to be of a normal kind. The 
distributions of the heat fluxes in the system and of the heat out-fluxes are of power-law shape. The basic 
properties of the model are in good qualitative agreement with experimental results. 
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I. INTRODUCTION 

Self-organized criticality (Soc) is a possible state of 
complex, spatially extended systems that are systemati- 
cally driven and that have mechanisms to develop local 
instabilities and to relax them. Soc is characterized by 
intermittent transport events that range from very small 
size up to system size, so called avalanches, by power- 
law distributions of variables characterizing the trans- 
port, and by power-spectra of the dissipated energy that 
are of power-law shape. Soc models are usually con- 
structed in the form of Cellular Automata (CA), i.e. by 
using a discrete grid and rules for the evolution of the 
system, with proto-type the sand-pile model of Ref. 0. 

In inquiries on confined plasmas and the related trans- 
port phenomena, evidence has been collected that the 
plasma might well be in the state of Soc, e.g. the fluc- 
tuations in density, potential, particle flux, and electron 
temperature show avalanche-like characteristics, such as 
intermittency and power-spectra of power-law shape, and 
the probability distribution of the particle flux displays 
a power-law (e.g. Ref s. [24 I27I . l26l and see the re- 
cent discussion in Ref. I33J ). Also fluid simulations have 
been found to exhibit Soc features, such as the oc- 
currence of avalanches and the appearance of frequency 
spectra of power-law shape (e.g. Refs. M, [32). A dis- 
crepancy between Soc models and experimental results 
had appeared concerning distributions of waiting-times 
between bursts in density-fluctuation and particle-flux 
times-series, the distributions determined from experi- 
mental time-series were found to be of power-law shape 22 , 



whereas the earlier Soc models, close in form to the 
model of Ref. 2s, yielded distributions of exponential 
shape. It was though realized later that spatial corre- 
lations in the loading process can lead to power-law dis- 
tributed waiting-time o 1 ' 29 ' 33 , and Ref. 23 found that, de- 
pending on the threshold applied in the burst detection, 
the waiting time distributions can turn from exponential 
to power-law shape, even in the case of the classical Soc 
model of Ref. Q- Waiting time statistics can therefore 
not be considered as an appropriate tool to test for Soc 
behavior in physical systems. 

Several Soc models in the form of CA have been sug- 
gested for fusion plasmas that address different aspects 
of turbulent transport and that are able to reproduce 
a number of observed phenomena, including the trans- 
port suppressing role of sheared poloidal flows^i, the 
scaling of transport characteristics with system size^, 
the power-law shape of waiting time distribution s 29 ' 30 ' 33 , 
the occurrence of anomalous, super-diffusive transport 
events 5 , the appearance of enhanced confinement, edge 
pedestals, edge localized modes, and the L- to H-mode 
transitio n 6 ' 10 ' 11 ' 31 , and the non-diffusive energy trans- 
port observed in off-axis heating experiments 20 . 

These Soc models are all CA models of the sand- 
pile type, and they basically are variants of the original 
Soc model of Ref. @ and its generalizations by Ref. LL6|, 
which also allows non-local relaxation of the instabili- 
ties, and Ref. 1 12L which introduces the running sand-pile 
model that is continuously loaded. The basic elements 
of the models thus are sand-grains, height and height- 
differences of sand-columns, with a usually vague identi- 



fication of these system variables with the physical vari- 
ables such as density, energy density, or temperature, and 
with an alike vague association of the system dynamics 
with actual physical processes. (A different approach is 
followed in Ref. 35;, where a partial differential equations 
is constructed that exhibits characteristics of Soc, and 
in which natural variables can be used.) 

Here, we introduce a new Soc model in the form of 
a CA for ion temperature gradient (ITG) mode driven 
turbulence in fusion plasma. A main characteristic of 
the model is that it is constructed in terms of the usual 
physical variables and that the temporal evolution of the 
CA, which necessarily is in the form of rules, mimics 
actual physical processes as they are considered to be 
active in the system. All elements of the model are thus 
consistently interpretable in the usual physical way, no 
sand-grain or sand-pile analogy is used. 

In the low /3 core plasma of a tokamak, where the mag- 
netic field dominates over pressure, turbulence is driven 
by two main electrostatic micro-instabilities, the ITG 
driven modes and the Trapped Electron Modes (TEM; 
e.g. Ref. 0). In principal, the two instabilities can co- 
exist, in the case though we consider here, where ion 
heating dominates, the ITG modes become unstable and 
actually are the dominant instability. In the strong tur- 
bulence regime then, the radial ion temperature profiles 
stay close to marginally stable, the gradients are stuck 
to their critical values, which is termed profile stiffness 
or profile consistency (e.g. Refs. I7I.I28L IM). This plasma 
behaviour is very reminiscent of Soc, and it is a main 
purpose of the application shown in this article to ade- 
quately model it. 

From the CA and Soc modeling point of view, the 
necessary conditions for a system to be able to reach 
the state of Soc are: (i) there must be a driving pro- 
cess which systematically increases the 'stress' in the 
system, (ii) a threshold dependent instability must be 
defined in some way, and (iii) a relaxation process must 
set on if somewhere the instability threshold is exceeded. 
All these processes should preferably act only in a lo- 
cal neighbourhood, and they must be formulated in the 
form of discrete evolution rules. ITG mode driven tur- 
bulence obviously meets the three prerequisites for SOC: 
the system is systematically driven, namely heated, the 
ITG mode instability is threshold dependent, and there 
is obviously a process that relaxes the instabilities, given 
the fact that ITG mode driven turbulence reaches a sat- 
urated stationary state (e.g. Ref. 0). 

Sect. II presents the model and explains how it is de- 
rived from usual physical processes in the plasma, Sect. 
Ill contains the results, and a discussion and the conclu- 
sions are given in Sect. IV and V, respectively. 



such as the tokamak, with the modeled domain being 
[R — sa, R + sa], with R the major radius and s < 1, so 
that only the core region is taken into account. 

The basic scalar grid variable at the grid sites Xt € 
[R — sa, R + sa], for i = 1, . . . , L, is the local ion tem- 
perature Ti = T(xi), considered as a positive real num- 
ber (i.e. assuming also non-integer values, in contrast 
to some sand-pile models). The grid is assumed to be 
equi-spaced, with grid-spacing Ax = 2sa/(L — 1). The 
grid-spacing is not considered 'infinitesimal', as in the so- 
lution of differential equations, but it is finite, of the size 
of some smallest scale of interest, here of the size of a 
typical ion Larmor radius. 



A. Instability criterion 

As is well known, see e.g. Refs. |T7] and [jjl, the ITG 
mode instability has a critical dependence on the normal- 
ized scale length Rj Lt, where Lt is the ion temperature- 
gradient scale-length, 1/Lt ■= |VT|/T, in the sense that 
ITG mode instabilities are triggered if the condition 
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(1) 



is fulfilled 

In order to use the instability criterion in the CA 
model, we have to determine the temperature gradient 
VT at the grid sites X{ . We thereto interpolate T^ locally 
in the neighbourhood (i — 1, i, i 4- 1) around the central 
site i with a second order polynomial, and we differenti- 
ate the polynomial at the point Xi, which yields 
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as an approximation to VT/T (and which is identical 
with the expression that a central difference scheme ap- 
proximation for VT would yield). After all, the grid site 
i is considered to be unstable if 



R\Si\> R/L cl . it 



(3) 



B. Redistribution rules 



We assume a normal diffusive process to be triggered 
around the grid site where an instability occurs, i.e. we 
assume an evolution according to a simple Fokker-Planck 
equation without advective term, 
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nk B d t T = -d x q, 



(4) 



II. THE MODEL 

We consider a one dimensional grid along the direction 
of the minor radius a in a toroidal confinement device, 



with ks Boltzmann's constant, n the number density, 
and q the classical heat flux, 



q {cl) (x) = - X nk B d x T, 
where \ i s t ne heat diffusivity. 



(5) 



1. General form of the redistribution rules 

In the CA, the local diffusion process is formulated in 
terms of a set of redistribution rules, which describe how, 
in the case of an instability, the temperature in the local 
neighborhood is redistributed, the latter being defined 
here as consisting of the central unstable site, say i, and 
its two nearest neighbours, i — 1 and i + 1. The general 
form of the redistribution rules can be written as 

T+ =T t +F (T i - 1 ,Ti,T i+1 ;6i 1 R/L ai! t), (6a) 
T^ 1 =T i±1 +F±{T i _ 1 ,T i ,T i+1 ;5 i ,R/L clit ), (6b) 

where the values after redistribution are denoted with a 
plus sign ( + ) as superscript, and Fq, F± are the changes 
in temperature, here written with their possible general 
functional dependence. In the derivation of specific redis- 
tribution rules, we will take into account actual physical 
processes and restrictions, and we will also make some 
assumptions. 



2. Energy conserving local diffusion 

We demand the redistribution rules to conserve the 
total temperature of the three involved grid-sites. From 
a physical point of view, the energy should be conserved, 
i.e. the product of the local temperature and the local 
density, so that we actually make here the assumption 
that the density is constant, which can partly be justified 
by the small spatial extent of the local neighbourhood, 
which is of the order of three ion Larmor radii, and by the 
fact that the density in tokamaks usually is quite close 
to uniform. It is to note that the local conservation of 
the grid variable in the CA is an important property for 
the system to be able to reach the Soc state. ^From Eq. 
([6]), the condition for local temperature conservation can 
be written as 



F_ + F Q + F 4 



0. 



(7) 



Assuming the density to be constant in a local neigh- 
bourhood, Eq. (|4]) and ([5]) reduce to a simple diffusion 
equation, 



1 



-d t T = x d xx T, 



(8) 



where we have also assumed that the diffusivity x is con- 
stant over a local neighbourhood. Diffusion according 
to Eq. ([8]) results, under normal conditions, in the local 
smoothing of the temperature profile, and if the bound- 
ary conditions are fixed around the local neighbourhood, 
a linear asymptotic profile will be reached (in the one 
dimensional case we treat here), where diffusion stops 
since V 2 T = d^T = 0. The relaxation process in the 
CA should thus describe such a local profile smoothing, 
expressed in terms of rules, and for simplicity we assume 
here the smoothed profile to equal the asymptotic one 
and thus to be linear. 



The assumption for the relaxed profile to be linear to- 
gether with the assumption of temperature conservation 
already determine the temperature T 4 after relaxation 
at the central site i of an unstable region. Since the re- 
laxed profile is linear, T± equals the mean temperature 
in the neighbourhood, 
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C7£i+T+), 
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into which we insert Eq. (j6 
T i + F = ±(T i _ 1 

which we rearrange 

2F - F_ -F + = 



F_+T l+1 +F+), (10) 



-2T % + T,_! + T t 



i+l- 



(11) 
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Now we note that 2F - F_ - F + = 3F - [F 
F+], the sum in the square-bracket is though zero due to 
energy conservation [Eq. 0], and Eq. (jlll) determines 
the increment at the central site 
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For convenience, we define 



n := Ti - - (T,_i 



■Ti. 



(12) 
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In the following, the increments F± of the nearest 
neighbours will be derived from a condition on the lo- 
cal heat fluxes. 



3. Heat fluxes 

In the relaxation events described by Eq. ([6]), heat is 
transported in the local neighbourhood, and we quantify 
this heat transport by determining the associated heat- 
fluxes. The starting point is the definition of the dynamic 
heat flux, 
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W(x) = -n{x)k B T{x){v(x)), 



(14) 



where (v(x)) is a local average velocity with which heat 
is flowing. 

We consider a local neighbourhood i — 1, i, i-\- 1 around 
the grid site i, which is unstable. The mean location Si of 
the total heat content within the neighbourhood can be 
defined by using the temperature values Tj as weights, 



i+l 

E 3 AxT J 

E^ 

j=i— 1 



(15) 



and where we have made use of the assumption of con- 
stant density. Generally, s, will not coincide with a grid- 
point. During the relaxation event, the local heat content 
will move to a new mean location sf, so that the heat 
displacement is s~l — Sj. This displacement takes place in 
one time-step of duration At, so that the velocity with 
which heat is transported is 



m := (s+ - Si)/ At. 



(16) 



According to Eq. (fT4|). we can now define the heat flux 
of a relaxation event as 



1 i+1 

qt := -nk B Vi ^ T U 



(17) 



j=i-l 



or, with Eqs. flTBl and (TIB"]) and using the conservation of 
heat in relaxation events, J2)=i-i Tj = XiL-i F + , 
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l —nk B J2 j&Bpf-Tj). (18) 
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j=i-l 



Inserting the redistribution rules in the form of Eqs. 
into Eq. ([IS) yields 
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Q * = 2At nkB [{l ~ 1)F - + tF ° + {t + 1)F+] ' (19) 

where the sum in the square brackets equals £(F_ + Fq + 
F + )+F + -F-, which, with Eq. Q, reduces to F + -F- 7 
so that the heat flux finally writes as 



c l, 



(dy) 



Ax 
2At 



nk B [F+ - F_] 



(20) 



Local diffusion of normal nature 



writes as 

Ax , , , 1 

—-nk B \F+ - F] = — 

2At l + J 2Ax 



Xnk B (T i+1 -Ti-i), (23) 



in which, due to energy conservation [Eq. ([7)1], we can 
replace — F_ by F + F + , 



[2F+ + F ] = 



At 

(A^y 



- X (T i+1 -!■_!), (24) 



where Fo is given by Eq. ([12|). so that the final form of 
the increment F+ is 



F +=-l( F ° + T&r« T ™- T ". 



(25) 



With Eq. (J7J, the increment F_ follows as _F_ = — F — F + 
from Eq. (|2"5]). 
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F- = -- ( F - 7X — x ( Ti+1 - Ti-i; 



(26) 



and, for convenience, we define the normalized heat dif- 
fusivity a as 



At 

(A^y 2 



(27) 



We note that we do have physical units for Ax, which 
equals the ion Larmor radius, At though is in arbitrary 
units, as explained in Sect. Ill El below, so that also a 
necessarily is in arbitrary units. 



In order the local diffusion process to be of normal 
nature, we demand that heat can only be transported 
from a hotter to a colder site, since else we would intro- 
duce anomalous diffusive effects at the local scale. To 
ensure this condition, we demand that the dynamic heat 
flux q\ , Eq. (f2T3|). equals the classical heat flux qf of 
Fourier's law, Eq. ([5]). In this way, heat is transported al- 
ways in the direction of the driving gradient ('downhill'). 

We discretize the gradient in Eq. ([5]) with the same 
method as used in the derivation of the instability crite- 
rion in Sect. Ill At 
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-xnftB(2}+i-Ii_i). (21) 



With Eqs. (|20J) and (TSIJ, the condition 



(a) Unstable state. 



(b) Relaxed state. 



FIG. 1: (Color online.) Sketch of the relaxation process. 

The red hatched region in the unstable configuration 

(a) marks the excess amount of heat to be moved, and 

the blue hatched region in the relaxed state (b) marks 

the amount of heat received in the relaxation event. 

The green dashed line connects the temperatures at the 

grid sites and represents the local temperature profile. 



-,(dy) 



M) 



(22) 



5. Summary of the relaxation rules 

After all, the relaxation rules are (see Eqs. (|6j), (fT2j). 

(USD, tfHD, dMD, 423) 



T+ = T t - - Ti , (28a) 

3^i = r«±i + ~r, T ~(T i+1 - Ti-i), (28b) 

with free parameter the normalized heat difFusivity a > 
0, and where the increments Fq and F± in Eq. ([6]) are 
given as 



Ff) — ——Ti 



1 -T- ^ (T 



Ti_i). 



(29a) 
(29b) 



A sketch of the relaxation process is shown in Fig. [TJ 

The main effect of the local diffusion process in Eq. 
(|4|) should be to remove the cause of the diffusion, i.e. 
instabilities should be relaxed, and we numerically find 
that the redistribution rules we derived indeed relax the 
instabilities, i.e. it holds that R\5^\ < R/L clit , as long as 
the normalized difFusivity is positive and not too small, 
a > 0.01. 

Last, we note that the redistribution rules in Eq. (l28l) . 
without though the parameter a (i.e . formally for a = 0), 
were used in Refs. [la, LLSl and LL4J in the astrophysical 
context of Solar flares (in combination with an instability 
criterion different from Eq. (|3J|). 



C. Boundary conditions 

We assume a constant value T\, of the pedestal tem- 
perature outside the domain of the cellular automaton. 
Tfc is used in the calculation of the instability criterion, 
Eq. ([3]), and in the redistribution rules, Eq. (|28|) . at the 
edge of the system, formally entering the equations as 
To = Tl+i = TJ,. 

The redistribution rules together with the initial con- 
ditions (see below) ensure that there is no heat in-flow 
from the boundaries, since always Tb < T\ and T& < Tl, 
the temperature is always lower outside than inside the 
grid. 



(see the next section for a definition of the time-step), 
and (ii) we consider loading only in stable configurations, 
where there are no instabilities present in the system. 

Moreover, we consider different spatial heating pat- 
terns: global heating everywhere in the system, and heat- 
ing only in localized regions. The latter case includes 
heating only in the central region, reminiscent of Ohmic 
heating, heating only in a region off-axis, as in the case 
of rf heating, and a combination of central and off-axis 
heating, representing Ohmic and off-axis heating applied 
simultaneously. 

In either case of heating pattern, a grid site is chosen 
at random from the region that undergoes heating, and 
the temperature at this site is increased by a constant, 
predefined amount ASj . The value of the parameter ASj 
is chosen small enough to avoid over-driving the system, 
as usual for Soc models. 



E. Description of the algorithm 

The time evolution of the model is as follows: (i) In 
a first step, the grid is scanned, and a list of the un- 
stable sites is created, (ii) If there are unstable sites, 
the instabilities are relaxed during a second scan of the 
grid. Prior to redistribution, each site is checked again 
for stability and it is relaxed only if it is still found to be 
unstable. The reason for checking again is the following: 
It may happen that two neighbouring sites are simulta- 
neously unstable, the redistribution rules are then first 
applied to just one of them, whose relaxation may cause 
the other site not to be unstable anymore, and it would be 
physically unmotivated to redistribute its temperature. 

If the system is loaded only in the absence of instabil- 
ities, the steps (i) and (ii) are repeated until there are 
no instabilities in the system anymore, and only then the 
system is loaded again. Each completion of steps (i) and 
(ii) is considered a time-step, whose duration we con- 
sider to be At = 1 in arbitrary units. The repetition of 
steps (i) and (ii) , from the first appearance of an insta- 
bility until the system has reached an everywhere stable 
state, is termed an avalanche, as usual in the context of 
Soc models. Continuous loading of the system means 
that after the completion of a time-step the system is 
also loaded, independent of whether there are instabili- 
ties present in the system or not. In this case, there is 
no practical definition of avalanches possible. 



D. Driving process 

The system undergoes a driving (heating) process, in 
which the temperature is increased, simulating the effect 
of heat injection. The loading process as such is not 
further specified, and it is understood to represent Ohmic 
heating as well as external heating. 

Concerning the temporal characteristics of the loading 
process, we apply two variants, (i) we apply continuous 
loading, where the system is heated in every time-step 



III. RESULTS 



A. Parameters 



The parameters of the model are chosen to coincide 
with those of the Joint European Torus^ 7 (JET) exper- 
iment. We thus assume a minor radius a — 1.25 m and 
a major radius R = 2.95 m, and since the ITG instabil- 
ities are considered to dominate turbulence only in the 



core, we restrict the spatial domain used in the simula- 
tions to [R — 0.8a , R + 0.8a]. The grid-size Ax is of the 
same order as the ion Larmor radius pi, so that, with 
Pi ~ 0.49 cm under typical conditions for JET, we use a 
grid of L — 401 sites with grid-size Aa; = 0.5 cm. 

According to Ref. Qj| the threshold R/L CTit of the ITG 
mode instability can assume values in the range [3.5,5]. 
In the following, and unless stated otherwise, the insta- 
bility threshold is set to R/L crit = 4, the normalized dif- 
fusivity to a = 0.5, the heating increment to ASj = 0.5, 
and we apply a constant value of the pedestal tempera- 
ture Tb = 500 outside the domain of the model as bound- 
ary condition. For the initial configuration we let the 
temperature everywhere equal T\,. 

The basic free parameters of the model, which mainly 
determine the results and which can be used to adjust the 
model to experimental data, are the threshold i?/L cr it, 
the pedestal temperature at the boundaries 7],, and the 
normalized diffusivity a. 
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FIG. 2: (Color online.) Spatial mean of (a) the 
temperature T and of (b) the absolute normalized scale 
length \R/Lt\ as a function of time, for heating only in 
stable states and for continuous heating, and with two 
different heating intensities. The red horizontal dashed 
line in (b) shows the instability threshold R/L cri t, while 
the arrows mark the transition of the system to the SOC 

state. 
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FIG. 3: (Color online.) Temperature profiles during the 

SOC state for (a) different values of R/L crlt and T b , and 

(b) for different values of a, all with loading only in 

stable configurations. 



B. The Soc state 

Fig. [2] shows the spatial mean of the temperature T(t) 
and of the absolute local normalized scale length \R/Lt\ 
for typical simulations with global heating, one with load- 
ing only in stable configurations and one with continuous 
loading, and for two different heating intensities. There is 
a transient phase, in which both mean values are growing, 
and after which there is a turn-over to a dynamic equi- 
librium state, where the mean values fluctuate around an 
asymptotic value. Once the asymptotic state is reached, 
avalanches of widely varying sizes start to appear, from 
very small ones, bounded to a local region, to avalanches 
that sweep through t he entir e system. The reaching of an 
asymptotic value of \R/Lt\ and T(t), together with the 
appearance of avalanches of all sizes, are indicative for 
the system to have reached the Soc state. The duration 
of the transient phase depends on the intensity of the 
driving process. It is to note that, in the case of intense 
continuous loading, the mean normalized scale length can 
also be above critical, whereas in the case of loading only 
in stable configurations it is always below critical. 

Temperature profiles in the Soc state are shown in 
Fig. Efa) for different values of the threshold i?/L C rit 



and pedestal temperature T& (here and in the following, 
the right halves of the symmetric profiles are shown). 
In all cases, the profiles are of exponential shape, and 
the threshold, together with the value T{, of the pedestal 
temperature at the boundaries, determine the maximum 
value that is reached in the center: the central tem- 
perature increases with increasing threshold as well as 
with increasing pedestal temperature, the latter more- 
over causing a shift of the entire temperature profile to- 
wards higher values. 

During the Soc state, the dynamically evolving tem- 
perature profiles stay very close to the characteristic ex- 
ponential shape, even the distortions caused by large 
avalanches would hardly be visible in the scales of Fig. |3] 
This implies that the total energy content of the sys- 
tem is, within small fluctuations, constant, heat is ejected 
from the system at the boundaries at the same statistical 
rate it is injected by heating. 

In Fig. EJb), we show the temperature profiles during 
Soc state for different values of the normalized diffusiv- 
ity a, keeping the pedestal temperature and the thresh- 
old fixed. The central peak temperature decreases with 
increasing diffusivity, a large a causes heat to be faster 
redistributed in the system and also to be more efficiently 
ejected at the edge. 



C. Central and off-axis heating 

In order to investigate the impact of off-axis heating 
on the temperature profile, we apply two spatial heat- 
ing patterns, one for pure central heating and one for 
simultaneous central heating and equally strong off-axis 
heating (the central region being defined as the grid- 
sites i € [151,251], and the off-axis region as the sites 
i G [351, 401]). In Fig.lH the temperature profiles yielded 
by the two simulations are shown during the Soc state. 
Also the profile from off-axis heating is clearly peaked at 
the center, there is no sign of the off-axis heat source visi- 
ble, and the difference between the two profiles is actually 
very small. This implies that the temperature profiles ex- 
hibit a very high degree of stiffness, they are basically not 
influenced by the applied spatial heating pattern. This 
also implies that there is strong anomalous transport, 
where heat is transported against the driving gradient 
('uphill'), despite the fact that diffusion at the local level 
is imposed to be of a normal kind. 



D. Continuous loading 

In the cases of temperature profiles shown so-far, heat- 
ing was applied only when there were no instabilities 
present in the system. This is justified physically only if 
the ITG mode instabilities together with their relaxation 
have a much faster time-scale than the heating process. 
When the two time-scales are comparable, loading should 
take place more frequently. 



1600 K 



1400- 



Loading pattern 
-Central, continuous 

Central, only in stable states 
-Central and off-axis, only in stable states 



600- 



FIG. 4: (Color online.) Temperature profiles for central 

heating (black dash-dotted), and for simultaneous 

central and off-axis heating (red dashed), both with 

heating only in stable states, and for central, continuous 

heating (blue, solid). 



In Fig.[4j we also show a temperature profile for contin- 
uous heating in the central region, here heating in every 
time-step. The profile is also of exponential shape, and 
it is similar to the one with heating only in stable states, 
the temperature values are though slightly larger, be- 
cause there is more heat injected into the system in this 
case, and the system has less time to eject heat at the 
edges (see also Fig. [2]). In the case of continuous heating, 
the value of the central peak temperature depends also 
on the heating increments ASj , it increases with increas- 
ing ASj, whereas in the case of loading only in stable 
configurations it is almost independent of ASj . 



Distribution of heat out-fluxes 



The local heat flux per relaxation event is given by 

q t in Eq. (f20|) . so that the total heat flux in a time- 
step (see Sect. Ill Ej) is 



Qt :=^gi(t), 



(30) 



and the total internal heat flux of an avalanche qi n t fol- 
lows as 



lint 



= E Q* 

t £ avalanche 



(31) 



Correspondingly, the total heat out-flux per avalanche 
at the right edge qo.8 is given by setting i = L in Eq. 
(I2U1) (i.e. considering the radial position x = 0.8a) and 
summing over the avalanche, 



9o.8 := 



E «*(*)• 



(32) 



t^avalanche 
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FIG. 5: (Color online.) Distributions (normalized histograms) of the total heat fluxes per avalanche qi n t [(a), (c), 

(e)] and total heat-outfluxes per avalanche go. 8 [(b), (d), (f)], for different loading intensities ASj [(a), (b)], for 

different thresholds R/L clit [(c), (d)], and for different diffusivities a [(e), (f)], together with power-law fits (dashed 

lines). (The distributions are shifted in the vertical direction for better visualization, and in each sub-figure 

error-bars (red) are shown for reference in one case.) 



For the case of global heating in stable configurations, 
Fig. [5] shows the normalized distributions of the inter- 
nal fluxes qi n t and the heat out-fluxes go. 8 from a large 
number of avalanches, and for different values of the pa- 
rameters ASj, i?/L cr jt, and a. All the distributions are 
of clear single or double power-law shapes in the inter- 
mediate range, extending over roughly 3 decades for qi n t 
and 1.5 decades for go. 8- The larger extent of the power- 



laws of qi n t reflects the fact that, from its definition, q^ n t 
has a larger dynamic range. The power-laws we find, to- 
gether with Fig. [3J are indicative for the system to be in 
the state of Soc. 

The indices a and /3 of the fitted power-laws {p(qint) °c 
lint an d p(Qo.s) c* Qqs j respectively) are summarized in 
Table |U (1) The power-law indices are basically indepen- 
dent of the threshold R/L cllt . (2) They depend on the 



TABLE I: Power-law indices a and (3 from the power-law fits to the distributions of the total internal fluxes qi n t and 

the total out-fluxes &o. s, respectively, in Fig. [5j 
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heating intensity ASj, the dependence is though very 
weak if the ASj are small, i.e. for low-level heating, the 
latter being defined in the sense that the system, when 
in Soc state and just having become stable, must on 
the average be heated with much more than just one 
heat increment ASj in order to get unstable again. (3) 
The power-law indices clearly depend on the diffusivity, 
they decrease and the distributions become flatter with 
increasing <r, a high diffusivity at the local level thus 
favours large heat-fluxes at the macroscopic scale, as one 
also would intuitively expect. 



IV. DISCUSSION 

A. Marginally stable profiles 

^From Eq. ([I]), we can analytically determine 
the marginally stable configuration, defined through 
\VT(x)\/T(x) = 1/L cr jt, which has the general solution 
T{x) = ce _ ' x '' <="*, and where c depends on the bound- 
ary conditions we apply, T(x = £) = Tb (with £ := 0.8a), 
from where c = Tf,e £ / Lcrit , so that 



T(x) = T h e^N^/ L c 



(33) 



(An exponential, marginally stable profile is also derived 
in Ref. :8j, in the frame of a critical gradient transport 
model.) The dynamic profile of the model in the Soc 
state follows the shape of the marginally stable profile, 
with a smaller decay constant though, i.e. reaching lower 
temperature and thus being sub-marginal when the sys- 
tem is loaded only during stable states, and with a larger 
than critical decay constant in the case of intense enough 
continuous loading (see Fig. [2J . The maximum temper- 
ature in the centre is T(x = 0) = Ti,e^' L crit according 
to Eq. (I33|) . and in the numerical simulations we indeed 
find the peak temperature to be proportional to T), and 
1/^crit (see Fig.EJa)), and to be dependent on i such that 
larger systems reach larger temperatures in the center. 

It is worthwhile noting that Eq. ([3"3| does not give a 
complete qualitative description of the system dynamics, 
it does not contain any information on the local relax- 
ation process and therewith on the third free parameter, 
the normalized diffusivity a. The latter has though a 
decisive influence on e.g. the central peak temperature 
reached, as Fig. [3]shows, the peak temperature increases 



with decreasing cr, and the profile actually approaches 
the marginally stable profile for very small a. 



B. Comparison with experimental data 

Experimental ion temperature profiles in tokamaks are 
well known to be very stiff and of exponential shape in 
the inner core region (see e.g. Ref. Q). Examples of pub- 
lished, exponentially shaped ion temperature profiles in- 
clude different devices and confinement modes, e.g. the 
L-mode in Tokamak Fusion Test Reactor, TFTR 38 (Fig. 
1 in Refjp), the H-mode in ASDEX Upgrade^ (Fig. 2 
in Ref. i^7 the L-mode in ASDEX Upgrade (Fig. 2 in 
Ref. EL and the L-mode in DIII-D^ tokamak (Fig. 1(c) 
in Ref. Q), and the high degree of stiffness in experimen- 
tal ion temperature profiles is presented and analyzed 
in e.g. Ref. [jjl Moreover, Ref. |27| experimentally finds 
power-law distributions for particle fluxes (with power- 
law index 1), which can be related to heat fluxes under 
the assumption that the heat-flux is convective. 

We thus can conclude that generally a good qualitative 
agreement of the model with experiments can be achieved 
with respect to (i) the exponential profile shape (Sect. 
llTTBJ) . (ii) the high profile stiffness (Sect. IfflCj) . and to 
some degree there is also a qualitative agreement with 
respect to (iii) the power-law distribution of out-fluxes 
(Sect. IIIIE[) . Last, we note that, having used parame- 
ters of the JET experiment (concerning the system size, 
ion Larmor radius, and pedestal temperature Tb) and as- 
suming Tb to be in units of eV, (iv) we naturally find a 
dynamic range of the ion temperature that is comparable 
to the one experimentally seen at JET, see e.g. Fig. 2(e) 
in Ref. [l9| in comparison with our Fig. [3] For a quantita- 
tive comparison, we would need to find a way to calibrate 
the time and therewith the normalized diffusivity in the 
model. 



C. Remarks on the SOC modeling approach 

The purpose of the study presented here was to in- 
troduce a Soc model for the ion temperature dynamics 
in the core region of tokamaks, with the elements of the 
model being physically interpretable in a sound and con- 
sistent way. This implies that in the derivation of the 
instability criterion, relaxation rules, and loading pro- 
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cess, we had to be guided by and take into account the 
basic physical variables and processes that are active in 
the system (to the degree that they are known). 

The Soc model was constructed from four basic pieces 
of information: (i) the form of the instability crite- 
rion, Eq. (JT|), which was taken from ex peri ments on ITG 
mode driven turbulence (Fig. 1 in Ref. ll9T) . (ii) a simpli- 
fied Fokker Planck equation from general transport the- 
ory, Eqs. (QJ and jSJ, (iii) the definition of the classical 
(Fourier's law, Eq. ([5])) and dynamical (Eq. (ITU) ) heat 
flux, and (iv) the concept of Soc in its operational defi- 
nition (basically following Ref. 0, and also Ref. fl8l ). 

The operational definition of Soc actually expresses 
the necessary conditions for an extended system to reach 
the Soc state, which partly were already mentioned 
in the Introduction. An extended system disposed to 
Soc must (a) systematically be driven, (b) allow local, 
threshold-dependent instabilities, (c) posses a mechanism 
that relaxes the instabilities locally, and (d) the system 
variable considered must be conserved in the relaxation 
events. In the model construction, the concept of Soc 
played the guiding role, in that it determines the elements 
of the model and their form, it contains the formal pre- 
requisites for the final model to be a realization of the 
Soc concept. 

In the classical sand-pile model of Ref. 2 and in many of 
its descendants mentioned in the introduction, instabili- 
ties occur if the height difference between neighbouring 
sand-columns exceeds a threshold, where the height dif- 
ference can be interpreted as a very rough approximation 
to the gradient V/i of the local sand-column height h. In 
the tokamak core plasma though, instabilities occur when 
S7T/T exceeds a threshold (see Sect. lTTA")) . not VT. The 
core plasma thus belongs to a formally different class of 
unstable systems than the sand-pile, which implies that 
we cannot make use of the sand-pile paradigm in core 
plasma Soc modeling, as long as it is our explicite aim 
to incorporate the actual physics of the problem under 
consideration as close as possible into the Soc model. 

It is worthwhile noting that Soc models that belong to 
different instability classes exhibit different characteristic 
spatial profiles of the grid variable. The sand-pile mod- 
els, with an instability criterion formulated in terms of 
an approximation to the gradient V/i of the sand-column 
height, exhibit linear profiles, see Ref.@- In an astrophys- 
ical context, Refs. [lg, [l3|, and [lj use an approximation 
to the Laplacian, V 2 A, in the instability criterion, with 
A the vector potential of the magnetic field, and they 
find parabolic spatial profiles for the components of A. 
The instability criterion used here, an approximation to 
VT/T, leads to exponential spatial profiles, as shown in 
Sect. IIIII The reason for these characteristic differences 
is that, in all cases, the spatial profiles in Soc state follow 
the marginally stable profiles in shape, which are differ- 
ent for the different classes of instabilities. 

Last, we note that our model has a high computational 
efficiency, which is a great advantage over e.g. the gyro- 
kinetic approach, it allows for fast interpretation and 



analysis of experimental results, and, after some further 
development, it even might be of potential use for plasma 
control. For example, the temperature profiles can easily 
be modeled globally, in the entire core region, heating 
patterns of any kind can easily be applied and explored, 
and the statistics of the heat fluxes can be determined 
very accurately, since the model can be monitored over 
very long times. 



V. CONCLUSION 

We introduced a cellular automaton model that im- 
plements the basic physics of ITG driven turbulence, as 
relevant for the core of fusion plasmas in tokamak de- 
vices, in L-mode as well as in H-mode. The model is 
formulated in terms of evolution rules, which makes it 
computationally very efficient. The rules implement two 
processes, heating by simple heat deposition, and local 
diffusion if a threshold in the normalized ion tempera- 
ture gradient R/ Lj< is exceeded, whereby it is ensured 
that energy is conserved and that heat is never trans- 
ported from a cooler to a hotter site. The model is for- 
mulated in terms of the usual physical variables, also in 
what the threshold is concerned, and the actual physical 
processes are directly mimicked in terms of rules. 

The model also represents an implementation of Soc, 
which it always reaches after an initial transient phase. 
In the Soc state then, the model yields symmetric ion 
temperature profiles of exponential shape. These pro- 
files exhibit very high stiffness, in that they basically are 
independent of the loading pattern applied (central and 
off-axis heating yield the same profile). This implies that 
there is anomalous heat transport ('uphill' heat trans- 
port, against the driving gradient) present in the system, 
despite the fact that diffusion at the local level is imposed 
to be of a normal kind. 

In a qualitative comparison of the model's basic prop- 
erties with experimental data, we find good agreement, 
at least for time instances where the experimental pro- 
files also exhibit an exponential shape. We thus can con- 
clude that the physical system we investigate, i.e. ITG 
driven turbulence, is qualitatively compatible with the 
Soc state of our model. 
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